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Rotating spiral and scroll waves (vortices) are investigated in the FitzHugh-Nagumo model of excitable media. The 
focus is on a parameter region in which there exists bistability between alternative stable vortices with distinct periods. 
Response Functions are used to predict the filament tension of the alternative scrolls and it is shown that the slow- 
period scroll has negative filament tension, while the filament tension of the fast-period scroll changes sign within 
a hysteresis loop. The predictions are confirmed by direct simulations. Further investigations show that the slow- 
period scrolls display features similar to delayed after-depolarisation (DAD) and tend to develop into turbulence similar 
to Ventricular Fibrillation (VF). Scrolls with positive filament tension collapse or stabilize, similar to monomorphic 
Ventricular Tachycardia (VT). Perturbations, such as boundary interaction or shock stimulus, can convert the vortex 
with negative filament tension into the vortex with positive filament tension. This may correspond to transition from 
VF to VT unrelated to pinning. 
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Orderly contraction of the heart is essential to pump blood 
efficiently. This order is imposed by electrical excitation 
waves propagating throughout the heart muscle. Abnor- 
malities in the propagation of excitation waves, known as 
arrhythmias, are responsible for many cardiac patholo- 
gies. Particularly dangerous are re-entrant arrhythmias, 
where excitation waves circulate along closed pathways. 
Re-entrant excitation "vortices", that are not attached to 
anatomical features, are spiral waves (in two dimensions) 
scroll waves (in three dimensions). As a rule, the rotation 
frequency and shape of a spiral wave is uniquely deter- 
mined by local properties of the tissue. We investigate a re- 
gion of parameters in a simplified mathematical model of 
excitation, in which two classes of vortices can exist. They 
differ significantly in their frequency and shape. Which 
of the two is realized, depends on the initial conditions. 
The defining feature of the model, responsible for this di- 
chotomy, is well known in heart electrophysiology. It is 
called delayed after-depolarization and is known to have 
a pro-arrhythmic effect. We focus especially on the three- 
dimensional scroll waves. Based on the asymptotic theory 
utlizing Response Functions, we predict a change of sign of 
the filament tension of the scrolls on one of the branches. 
Negative tension is known to promote "scroll wave turbu- 
lence". We investigate, by numerical simulations, specifics 
of this phenomenon in presence of alternative scrolls. We 
describe conversion of one type of vortex to the other un- 
der various perturbations, including external shocks such 
as are used for defibrillation, interactions with boundaries, 
and curvature of the scroll filaments. Finally, we discuss 
implications of our findings for cardiac arrhythmias. 



I. INTRODUCTION 

Spiral waves in two-dimensions, and scroll waves in 
three-dimension, are r egim es of self-organization observed 
in physical^ chemicar'^, and biological^ dissipative sys- 
tems, where wave propagation is supported by a source of en- 
ergy stored in the medium. The interest in the dynamics of 
these waves has significantly broadened over the years as de- 
velopments in experimental techniques have permitted them 
to be observed and studied in an ever increasing number of di- 
verse systems^^ However, the occurrence of these waves in 
excitable media, and cardiac tissue in particular, has been and 
continues to be one of the main motivating factors for their 
study. 

In most situations, a spiral wave in excitable media rotates 
with a period determined uniquely by the medium. However, 
there are cases in which there is bistability between alternative 
waves in the same medium. It is precisely this case which is 
of interest in this paper. Our particular concerns are (1) the 
transitions between the alternative solutions, that is, how can 
transition from one type of solution to the other be effected, 
(2) the differences in the dynamics exhibited by alternative 
scroll waves in the 3D setting, and finally (3) the relationship 
these phenomena may have to cardiac electrophysiology. 

Our study is based on the FitzHugh-Nagumo (FHN) model, 
a standard two-component reaction-diffusion model capturing 
the essential features of excitable media such as cardiac tissue. 
The model is given by 

dtu = f{u,v)^\/^u, (1) 
dtv = g{u,v), (2) 



where u and v are the dependent variables, corresponding 
roughly to membrane voltage and ionic channels, respectively. 
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FIG. 1. (Color online) A pair of alternative vortices for the same 
model parameters: a — 0.3, /3 = 0.72, 7 = 0.5. (Left) slow 
solution and (right) fast solution. (a,b) spiral waves, shown is the 
i^-field, and tip trace; (c,d) action potentials; (e,f) phase portraits, 
enlarged is the area around the fixed point. 

The kinetic terms in the model are 

f{u^ v) = ar^{u — /?> — v), 
g{u, v) = a{u-\- f3 — jv). 

with parameters a, (3, and 7. Space units are set such that the 
diffusion coefficient is 1 . 

Typical solutions to the FHN model include rotating spiral 
waves in 2D and rotating scroll waves in 3D. These waves, 
whether in 2D or 3D, are commonly referred to vortices, al- 
though they are unrelated to fluid vorticity. As already noted, 
in most situations vortices rotate with a period determined 
uniquely by the model parameters. However, it has been 
shown by Winfree^^ that for certain model parameters there 
exist alternative stable spirals with distinct periods. A pair 
of such stable spiral solutions in the FHN model is shown in 
Fig.[T] We refer to the two solutions as slow and fast. The slow 
spiral has a longer temporal period, larger spatial wavelength, 
and larger core radius than the fast one. The longer period of a 
slow vortex is due to the small extra loop near the fixed point 
in the u-v phase portrait. Fig. ^0), which manifests itself as 
extra maximum in the tail of the action potential, Fig.[TJc): the 
feature known as delayed after-depolarisation (DAD) in car- 
diac electrophysiolog>^^ The fast, shorter period vortices do 
not show DADs in their action potentials, Figs.[TJd) and[TJf). 

The shaded region in Fig. [2] shows the portion of parame- 
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FIG. 2. (Color online) Portion of the parameter space of the 
FitzHugh-Nagumo model containing alternative vortices. 7 is fixed 
at 0.5. Within the shaded cusp-shapped region both fast and slow 
vortices exist. The thin vertical line at a = 0.3 shows the parameter 
cut considered throughout the paper. The dashed curve marks a line 
of zero filament tension. To the lower left of this curve fast vortices 
have positive filament tension. 



ter space, for fixed 7 = 0.5, where alternative vortices exist. 
This region has been computed as part of the present study. 
Fast spirals exist to the lower left while slow spirals exist to 
the upper right. Within the cusp-shaped region both types of 
vortices exist with bistability between them. The boundaries 
of the bistable region are fold singularies (limit points) meet- 
ing at a cusp, (a, /3) ^ (0.27, 0.8). Above the cusp, fast and 
slow vortices are connected continuously. The thin vertical 
line shows a representative one-parameter cut, a = 0.3, that 
will be the focus of much of our study. The dashed line relates 
to a 3D phenomenon which we now address. 

In 3D, scroll waves are organized about filaments and the 
possible behavior is richer than in 2D. Filaments are not, in 
general, fixed in space but instead undergo motion, typically 
on a slow timescale relative to the rotation period. Hence, in 
addition to whatever dynamics 2D spiral waves might have, 
scroll waves exhibit additional dynamics associated with fila- 
ment motion^^ Working in Frenet coordinates, the motion 
may be conveniently expressed in terms of the velocities V/v 
and Vb in the normal and binormal directions, respectively, 
at each point along the filament. Motion along the tangential 
direction is of no physical significance and is equivalent to 
reparametrization of the filament. 

First semi-phenomenologicall>^ and later using asymp- 
toticJ^^EH^ the equations of filament motion have been ob- 
tained. At lowest order these are 

Vat = 62/^ + . . . , VB=csf^^..., (3) 

where k, is the filament curvature. The coefficients 62 and C3 
depend on properties of the medium and 2D spiral solutions. 
The coefficient 62 is called the filament tension. To understand 
this, consider a circular filament, i.e. a scroll ring. For posi- 
tive 62 a scroll ring will contract somewhat as if the ring were 
an elastic ring under tension. For negative 62, the ring will 
expand, as though under negative tension. More generally in 
the negative-tension case, filaments will increase in length and 
eventually evolve into full scale autowave turbulence^^. The 
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FIG. 3. (Color online) Alternative vortices along the parameter cut 
a — 0.30, 7 = 0.5 shown in Fig. [2] Red squares are fast spirals 
and blue circle are slow spirals, (a) Hysteresis loop in the rotation 
velocity u. (b) Prediction for filament tension 62. (c) Prediction for 
the binormal drift coefficient C3. 



dashed curve in Fig. |2] indicates a change of sign of filament 
tension discussed at length later in the paper. The coefficient 
C3, which we shall call the binormal drift coefficient describes 
the drift of a scroll ring perpendicular to the plane of the ring, 
or more generally, the velocity component orthogonal to the 
local plane of the filament. 

The remainder of the paper is then devoted to understanding 
the alternative vortices in 2D and 3D. We first consider rele- 
vant asymptotic theory and discuss how Response Functions 
can be used to predict the filament tension of the alternative 
vortices. We then investigate the dynamics of the alternative 
vortices in 2D and 3D through direct numerical simulations. 
In addition to confirming the prediction, we study how the 
alternative vortices can be converted one into another in a va- 
riety of situations. 

II. ASYMPTOTICAL PREDICTIONS 

A. Filament tension and drift 

At leading asymptotic order, the motion of scroll filaments 
is determined by the two coefficients 62 and C3 appearing in 
Eqs. ([3]). Once these coefficients are known, the most impor- 
tant properties of filament dynamics, especially the sign of the 
coefficient 62 determining the sign of the filament tension, are 
determined. 

The coefficients for filament motion in Eqs. ([3]) are given 
by the following simple formulsP*^ 
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FIG. 4. (color online) Spiral waves, shown is the i^-field (below, flat), 
and their translational response functions jwi^'*! (above, elevated 
surface) for (a) the slow solution, and (b) the fast solution; /3 = 0.72. 



stone mode, W^^^ is the corresponding Response Function 
(RF), and D is a diffusion matrix. We now elaborate, although 
we refer th e reader to other publications for most technical de- 

First consider a rotating spiral wave solution of Eqs. ([TJ-Q- 
It is convenient to work in a frame of reference co-rotating 
with the spiral, at angular velocity uo. Denote the steady spiral 
seen in this co-rotating frame by U = (li, v)^ . 

Now consider the linear stability of such a spiral. Due 
to symmetries of the system, there will be three eigenvalues 
given by A = mcj, with corresponding eigenvectors or Gold- 
stone modes V'^^^ where n = — 1, 0, 1. The n = eigenvalue 
and Goldstone mode is related to rotational symmetry and the 
complex pair of eigenvalues and Goldstone modes n = ±1 is 
related to translational symmetry. 

In addition to the linear stability problem there is the asso- 
ciated adjoint problem. The adjoints eigenmodes correspond- 
ing to eigenvalues —inuj are the RFs and are denoted W^^^. 
Thus one sees that to obtain the coefficients for filament mo- 
tion in Eq. ^ we require the translational Goldstone mode 
V*^^^ and corresponding Response Function W*^^^. D is just 
the matrix of diffusion coefficients appearing in the reaction- 
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where angle brackets denote an inner product, V*^^^ is a Gold- 



diffusion equations, so for Eqs. ( 1 

angle brackets in Eq. ^ signify a straightforward integration 
over space of the Hermitian product of the vector fields. We 
note that similar methods employing RFs can be used to ob- 
tain the drift velocities of vortices in response to perturbations 
to Eqs. Q-llPt^. 

The prediction of filament tension, and more generally drift 
velocities, relies on two conditions: firstly, the localization 
of the RFs W*^^^ in the vicinity of the core of the spiral and 
secondly, the ability to compute the RFs efficiently and ac- 
curately. Existence of localized RFs has been demonstrated 
for a broad range of the models' parameters in several mod- 
gj^36 42 ^ robust method to compute the RFs with good ac- 
curacy for any model of excitable medium with differentiable 
kinetics has been developed in Ref. 41 , extending eariler sta- 
bility methods"^^. Using these methods, we have computed 
the steady spiral U, its rotational velocity uj together with the 
Goldstone modes V^"^^ and the response functions W^"^^ for 
model ([l])-([2]) on a disk of radius Pmax = 25 with Nq = 64 
grid points in the angular direction and Np = 251 grid points 
in the radial direction, see Fig.|4] 
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We focus now on the parameter path indicated in Fig. [2] 
For a = 0.3 and 7 = 0.5 the hysteresis loop in the spirals' 
rotational velocity uj has been obtained by continuation in pa- 
rameter (3. As is seen in Fig.[3ja), for these parameters, alter- 
native stable spiral wave solutions with distinct u exist in the 
range A < /3 < where /3e ^ 0.708, f3u ^ 0.740. 

We have also computed the filament tension 62 and the drift 
coefficient C3 along this parameter path and the results are 
shown in Figs.|3jb) andjSjc). The most significant feature of 
these results is the sign change of the filament tension within 
the hysteresis loop. Below (3 = (3^ ^ 0.722, the alterna- 
tive vortices have filament tension of opposite signs: the fast 
vortex has positive tension while the slow one has negative 
tension. At f3 = the filament tension of the fast vortex 
changes sign so that in the parameter range (3^ < (3 < Pu, the 
alternative vortices both have negative filament tension. 

The drift coefficient C3 computed in this parameter cut is of 
fixed sign (negative). The graphs of this drift coefficient of the 
alternative vortices cross Sit (3 ^ 0.726. 

B. Scroll rings and electrophoretic drift 

There is a connection between filament motion of scroll 
rings and the drift of spiral waves in response to applied elec- 
tric fields (electrophoretic drift). Both things will appear in 
the subsequent numerical studies and it is appropriate to sum- 
marize the issues here. 

A scroll ring not only has a circular filament, but the entire 
solution has an axial symmetry. It is convenient to study such 
a structure in cylindrical coordinates (r, 6>, z), where by sym- 
metry, the solution is independent of 0. Hence the diffusion 
term in Eq. ([T]) becomes 

In this way scroll rings may be studied using 2D numerical 
computations as long as one is not interested in any symmetry- 
breaking effects. Moreover, if one is interested in scroll rings 
with small curvature hz, corresponding to rings of large fila- 
ment radius R = 1/k,, one may use the approximation 

For small-curvature rings this is an exceedingly accurate ap- 
proximation. In fact it is equivalent to considering only the 
lowest-order curvature contributions as in Eqs. ([3]). We shall 
use this approach to evaluate filament tension and drift. Note 
that the normal to the scroll ring points in the outward radial 
direction and so a positive normal velocity Vn > corre- 
sponds to decreasing r, and vice versa. 

Returning to the 2D Cartesian situation, if the excitable 
medium is a reaction-diffusion system in which the u reagent 
is electrically charged and the v is neutral, then the effect of 
applying an electric field can be modeled by including a term 

4, (6) 

OX 



on the right hand side of Eq. ([T]), where E is proportional to the 
applied electric field which is taken to be in the x direction. 

For small values of £, the effect of such a term will be drift 
of the spiral and one can calculate, using Response Functions, 
the drift velocit>^^ However, one can use the formal equiv- 
alence^ between the term £du/dx of Eq. ^ and the term 
ndu/dr in Eq. ^ to obtain the drift velocity for vortex rings 
from the results for the applied electric field. Note, that elec- 
trophoretic drift towards negative x corresponds to positive 
filament tension and towards positive x corresponds to nega- 
tive filament tension. 



III. TWO-DIMENSIONAL SIMULATIONS 

2D simulations have been perfo rmed with a suitably mod- 
ified version of ezspiralHSI. 

Unless specified other- 
wise, simulations have been performed using forward Eu- 
ler timestepping on a uniform Cartesian grid on square do- 
mains 40 X 40 s.u. with non-flux boundary conditions, nine- 
point approximation of the Laplacian, spatial discretization 
Ax = A?/ = = 1/3 and time step At 3/80. The 
tip of the spiral is defined as the intersections of isolines 
u{x^y) = and v{x^y) = v^, and the angle between \/u 
at the tip and x axis is taken as orientation of the tip. We 
use (u^^v^) = (0, 0) for the FitzHugh-Nagumo model. Initial 
conditions for the alternative spirals have been computed us- 
ing parameter continuation and then converted into the input 
format of EZSPIRAL. 



A. Alternative vortices: conversion by a shock 

First, to test if the alternative vortices can be converted from 
one into another in a controllable way, we add a constant value 
A to the fast variable u, uniformily in space, at a specified time 
instant T. 

u{x,y^t) ^ u{x^y^t) -\- at t = T (7) 

This may be considered as a very crude model of a defibrillat- 
ing shock. We use a box size of 50 x 50 s.u., and in selected 
simulations up to 150 x 150 s.u., with spatial discretization 
h = 0.1 and time step At = 2.25 x 10"^. 

The minimum shock amplitudes required to convert one 
alternative vortex into another are shown in Fig. |5] together 
with the corresponding "defibrillation thresholds" that are suf- 
ficient for the complete elimination of spiral- wave activity. 

It can be seen that, within the parameter range of the hys- 
teresis Pi < (3 < there are three distinct intervals of be- 
havior demarked by the values Pi ^ 0.722 and P2 ^ 0.736. 
In the region between Pi and Pi it is possible to convert a slow 
vortex into its faster counterpart. Transitions in the opposite 
direction have not been observed, even as we increase A to 
the defibrillation threshold, which eliminates spirals from the 
medium. In the middle interval, between Pi and P2, conver- 
sion both ways is possible. For P in the range between p2 and 
Pu, we only observe conversions from fast vortices into their 
slow counterparts, but not vice versa. 
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FIG. 5. (Color online) Minimum shock amplitudes A needed to 
convert one type of spiral into the other (lines with large symbols), in 
comparison with the corresponding defibrillation thresholds (dashed 
and dotted lines). In the legend, 'conv fast'/' conv slow' refers to the 
shocks applied to fast/slow spirals necessary for conversion. Like- 
wise, 'defib fast'/'defib slow' refers to the shocks applied to fast/slow 
spirals necessary for defibrillation. The thick solid lines above the 
'conv slow' curve indicate the region of parameters causing conver- 
sion with multiplication via front break-up. 



The value of /3i , which is the smallest (3 for which conver- 
sion is possible from fast to slow vortices with a single shock, 
is close to (3^ at which the filament tension of the fast vortices 
is zero. This is a purely empirical observation and we have 
not investigated how precisely it holds. 

We also observed conversion with multiplication from a 
slow vortex with a large core into multiple fast vortices with 
small cores, via break-up of some segments of excitation 
fronts. This is illustrated in Fig. [6] This conversion with 
multiplication occurs in a very small range of shock ampli- 
tudes just above the conversion amplitude; the numerical val- 
ues are listed in Table |T| Note that break-up of spiral waves by 
spatially uniform shocks has been previously observed, e.g. 
Ref . |46j in cases without bistability, but at shock magnitudes 
close to the defibrillation threshold. Here we observe it in 
a system with alternative spiral waves, at shock magnitudes 
significantly smaller than the defibrillation threshold and only 
slightly above the conversion threshold. We have not observed 
conversion with multiplication from a fast spiral into multiple 
slow spirals by any shock amplitude. 



B. Conversion due to interaction witli the boundary 

To study how interaction with a boundary may bring about 
conversion, we gently push the spiral towards a boundary us- 



ing electrophoretic drift discussed in Sec. |IIB[ Recall that neg- 
ative tension (62 < 0) corresponds to electrophoretic drift to 
the right and positive tension (62 > 0) corresponds to drift to 
the left. 

Representative simulations are shown in Fig.|7]at /3 = 0.71. 




(a) 



(b) 



(c) 



FIG. 6. (Color online) Shock conversion with multiplication at 
/3 = 0.73. (a) Slow vortex with the large core before the shock, (b) 
Front breaks immediately after a shock of amplitude A = 0.29. (c) 
Multiple fast spirals with small cores, approximately 50.3 t.u. after 
the shock. 
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TABLE I. The 
tudes producing 


minimum Amin and maximum Amax 
I multiple spirals from a slow vortex. 


shock ampli- 



Simulations are started from the slow, large-core spiral. The 
period of the spiral is T ^ 16.8654. The spiral drifts from 
its initial position to the right in agreement with the predicted 
negative filament tension, 62 < 0. It also moves up, as it 
should, since predicted drift coefficient, C3 is negative. When 
the spiral nears a boundary, its core radius decreases signifi- 
cantly and its period changes to T ^ 11.2964. The horizon- 
tal drift component changes sign [Fig.jTJb)], corresponding to 
switching to positive filament tension, and the spiral moves to 
the left rather than to the right. The binormal component of 
drift changes magnitude but not sign. 

At the top boundary, the spiral drifts to the left until it pins 
to the top left comer of the box [Fig. [tJc)]. If the Neumann 
boundary conditions (NBCs) at the top and bottom edges of 
the box are changed to periodic boundary conditions (PBCs), 
the fast spiral with a small core will continue to drift upwards 
along the left boundary indefinitely [Fig.|7jd-f)]. 
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(a) (b) (c) (d) (e) (f) 

FIG. 7. (Color online) Conversion due to do interaction with a boundary. Simulations are started with a slow spiral at /3 = 0.71. The spiral is 
induced to drift toward the right boundary by electrophoresis perturbation of amplitude 8 = 0.03. (a-c) Three snapshots of the evolution with 
Neumann boundaries conditions on all four sides, (d-f) Three snapshots of the evolution with Neumann boundaries conditions on right and 
left, and periodic boundary conditions on top and bottom. 



C. Conversion due to applied field 

Electrophoretic driving of a sufficiently high magnitude can 
itself cause conversion of vortices directly. Figure [5] shows 
selected simulations illustrating this phenomenon. FigurejSJa) 
shows that 8 = 0.035 is not strong enough to convert the slow 
spiral into the fast, and indeed it allows the spiral to continue 
to drift to the right boundary, at which the conversion happens. 
Figure [8jb) shows the simulation with the threshold amplitude 
£ = 0.0472, allowing the conversion to happen well before 
reaching the boundary. Finally, Fig. [Sj^c) shows a simulation 
with £ = 0.06 which is strong enough to convert the slow 
spiral into the fast one almost immediately. 

Figures [8jd-f) illustrate the combined effect of an elec- 
trophoretic driving and defibrillation- style shock. Shocks of 
various amplitudes are applied after 8 rotation periods of the 
slow spiral. If the amplitude is sufficiently high, conver- 
sion happens before the spiral reaches the boundary. In the 
presence of electrophoretic driving with £ = 0.03, we find 
that the minimal shock strength required for conversion is 
^conv ^ 0.1489 (the case shown in the figure), which is 
slightly smaller than that Aconv ^ 0.1789 required for con- 
version in the absence of applied field (£ = 0). 



D. Verfication of predictions and effects of discretization 



both 62 and C3 converge extremely well as the grid size goes to 
zero to the predicted values given by Eq. Q. The convergence 
is quadratic, i.e. linear in h'^, as is consistent with the second- 
order accuracy of the numerical simulation. 

The case shown in the right column of Fig. [9j where the 
predicted value of 62 is small, is more problematic. Again, C3 
converge extremely well, and with the expected form, to the 
predicted result. The convergence of the tension 62 is less sat- 
isfying, both in the form of the convergence and in the asymp- 
totic value. The exact reasons for this are outside the focus of 
the current study. The important point is the sign change of 62 
seen in the right column of Fig.[9]due to finite resolution. This 
means that a 3D numerical simulation performed at even high 
resolution, e.g. h = 0.2 will be qualitatively different from a 
fully resolved simulation at these parameter values, since the 
simulated filament tension will be negative whereas the fully 
resolved tension is positive. 

The conclusion is that the filament tension and binormal 
drift predicted from RFs and plotted in Fig. [3] are borne out by 
direct numerical simulations. Nevertheless, 3D simulations 
should only be performed for the model parameters such that 
the discretization does not qualitatively affect the simulations 
by artificially changing the sign of the filament tension. For 
the parameters we consider, this means to the left (P < 0.71) 
or the right (P > 0.73) sides of zero crossing of 62 in Fig. [3] 



Before proceeding to 3D simulations, we compared the pre- 
dictions for the filament tension 62 and binormal drift coef- 
ficient C3 obtained via RFs in Sec. |ll) with direct 2D com- 
puter simulations of axisymmetric scroll rings as discussed in 



Sec. II B From the motion of the axisymmetric filament in the 
radial (normal) and vertical (binormal) directions we obtain 
the coefficients 62 and C3. In 2D, it is possible to perform sim- 
ulations to very high spatial resolution with reasonable cost 
and hence we are also able to study systematically the effects 
of numerical resolution. 

Simulations of electrophoretic drift at £ = 0.01 have been 
conducted for the fast vortex over a variety of regular grids 
with spacings Ar = Az = h. Results are shown in Fig.|9]for 
two values of /3: /3 = 0.71, a typical case away from where 
62 crosses zero in Fig.|3j and (3 = 0.722, a case very near the 
zero crossing of 62- In the typical case, left column of Fig.[9j 



IV. THREE-DIMENSIONAL SIMULATIONS 

We have seen in Sec.[IIl|that it is possible convert between 
fast and slow vortices in 2D by a variety of mechanisms, such 
as interactions with boundaries and shocks. These effects will 
also exist in 3D, but in 3D the additional effects arising due 
to filament curvature and filament tension can be hugely im- 
portant. This is especially so in the region where, according 
to predictions of Sec.jllj the alternative vortices have opposite 
signs of filament tension. 

Simulations have been performed with a version of 
EZSCROLL47 48 modified for FHN kinetics. The simulations 
use forward Euler timestepping on a uniform Cartesian grid on 
cuboid domains of varying size with non-flux boundary con- 
ditions, nineteen-point approximation of the Laplacian, space 
discretization step Ax = Ay = Az = h = 1/3 and time 



a) (b) (c) (d) (e) (f) 

FIG. 8. (Color online) Conversion due to allied fields. All simulations start from a slow spiral at /3 = 0.71. (a-c) Three simulations for field 
strengths (a) 8 = 0.035, (b) 8 = 0.0472, (c) 8 = 0.06. In (a) conversion occurs at the boundary, whereas for (b,c) conversion occurs due 
to sufficiently strong electrophoretic drift, (d-f) Three snapshots of a single simulation with 8 = 0.03. After 8 rotations a shock of strength 
A — 0.14889 is applied resulting in conversion. 
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FIG. 9. (Color online) Numerical convergence of 62 and C3 from 
numerical simulations of electrophoretic drift. The fast vortex is 
simulated over regular grid spacings Ar — Az — /i as indicated. 
Left column: /3 = 0.71, a typical case away from 62 = 0. Both 
coefficients converge quadratically in h to the values predicted by 
RFs (blue dashed lines). Right column: P = 0.722, a case close to 
62 = where filament tension is small. C3 converges quadratically 
while the convergence of 62 is less clear. Importantly, except at very 
high resolution the tension in the simulation is of the opposite sign 
from the asymptotic result. 



discretiation step At = 3/80. Numerical determination of the 
instant position of scroll filaments is technically challenging, 
and we use an easy substitute: the instant phase singularity 
lines, defined as the intersections of isolines u{x^y) = 
and v{x,y) = v^. Such singularity lines correspond to spi- 
ral tips in 2D, and like a spiral tip rotates around the current 
rotation center, so a singularity line rotates around the current 
filament. We also assume that inasmuch as the twist of scroll 
waves is insignificant, then the curvature of the singularity line 
is close to the curvature of the filament, so we take the for- 
mer as an approximation of the latter. As with 2D simulations 
in EZSPIRAL, we use {u^,v^) = (0,0) for the FitzHugh- 
Nagumo model, with some exceptions, described later. Initial 
conditions for alternative scrolls are obtained by first comput- 
ing alternative spirals on a polar grid, using methods described 
in Ref . ,41, then converting these solutions to Cartesian coor- 
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FIG. 10. (Color online) Evolution of a fast helical scroll with posi- 
tive filament tension at /3 = 0.71. The domain has Neumann bound- 
ary conditions, (a-f) The filament straightens up (the numbers un- 
der the snapshots show corresponding time in timeunits); (g) action 
potentials with blue crosses marking the times at which two of the 
snapshots are taken. Singularity lines are used as approximations to 
the scroll filament. 



dinates and "stacking" spirals on top of one other to generate 
3D scroll waves. 

It is more difficult in 3D than in 2D to use core size and vi- 
sual inspection to distinguish between fast and slow vortices. 
However, point records of the action potentials, u{t) at some 
fixed location, may clearly distinguish the two types of vor- 
tices. Recall Fig. [T] Slow vortices have extra maxima in the 
tails of their action potentials, the DADs, while fast vortices 
do not, no DADs. This method of distinguishing between vor- 
tices assumes that the whole solution consists of only slow or 
only fast vortices. A more sophisticated method free from this 
assumption is described later. 



A. Alternative vortices with opposite signs of filament 
tension 

The first case we consider is at /3 = 0.71, where the fast and 
slow vortices have opposite signs of filament tension. Simula- 
tions have been conducted in a box of size 50 x 50 x 50 s.u., 
with Neumann boundary conditions (NBCs) on all sides. Sim- 
ulations are started from alternative vortices with filaments in 
the form of helices constructed by layering 2D solutions such 
as to create a helical scroll with one turn from bottom to top. 
The radius of the helix is 2 s.u. We expect that the filament of 
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FIG. 11. (Color online) Evolution of a slow helical scroll with neg- 
ative filament tension ai (3 = 0.71. Same domain and parameters as 
for the fast vortex in Fig. [T0| (a-f) (Top) isosurfaces of the ix-field; 
(middle) filament only (as in Fig.[TO]). (g) action potentials with blue 
crosses marking the times at which the snapshots are taken. The helix 
initially expands as expected for negative tension, but then undergoes 
conversion and contracts. 



the fast vortex with positive filament tension will straighten up 
while slow vortex with negative filament tension will expand, 
possibly breaking up and developing into turbulence. 

Figure [T0| shows the evolution of the fast vortex with pos- 
itive filament tension. As expected, the filament straightens 
up, see Fig. [lot a-f). The record of corresponding action po- 
tential. Fig. [lO[g), shows that there is no conversion into the 
slow alternative vortex: no DADs appeared. 



Figure 1 1 shows the evolution of the slow vortex with neg- 
ative filament tension. First, the filament expands as expected. 
However, after five full rotations, the vortex spontaneously 
changes its period and converts into its fast counterpart with 
positive filament tension, and we see DADs disappearing at 
this moment in Fig. [TTJg). The fast, positive filament ten- 
sion vortex subsequently contracts, see Fig. [TTJd-f), and dis- 
appears. 



B. Conversion due to curvature 

To identify the cause for the spontaneous conversion in the 
simulations just discussed, we have repeated the simulation 
with the slow helical vortex in different domain sizes with pe- 
riodic boundary conditions (PBCs) on the top and bottom of 
the domain, and NBCs elsewhere. 

In this series of simulations, we start from a helical fila- 
ment with one full turn from bottom to top. We monitor the 
mean radius of the projected filament onto the base of the box. 
Initially the projection is a circle. The initial helical filament 
expands. As this expansion is inherently unstable, the pro- 
jections of the filament onto the base very soon cease to be 
circular. See Figs. T2p4 



Figure [T5[a) shows how the mean radius of a vortex fila- 
ment's projection changes every period in the boxes 50 x 50 x 
50 S.U., 50 X 50 X 100 s.u. and 70 x 70 x 50 s.u. The vertical 
lines show the time at which conversion happens in each box. 




(g) 



FIG. 12. (Color online) Conversion of a slow helical vortex into a 
fast one at /3 = 0.71. Domain size is 50 x 50 x 50 s.u. with PBC on 
top and bottom and NBC elsewhere, (a-f) evolution of the filament 
(top row) and projections of the filament onto the base of the box 
(second row) shown in red, with the barycenter shown as a blue cross 
and the mean radius shown as the green circle, (g) action potentials 
(i^-field against time t), blue crosses correspond to the times at which 
the snapshots in (a-f) are taken, and the period number is shown on 
the upper axis, (h) period T plotted against time t. 



as indicated by the morphology of the point records. It can be 
seen that in the boxes 50 x 50 x 50 s.u. and 70 x 70 x 50 s.u., 
both of height Lz — 50 s.u., conversion happens at approx- 
imately the same time, regardless of the width of the box. 
Whereas in the box 50 x 50 x 100 s.u. with the height lOOs.u. 
conversion happens much later than in the other two cases. 
Comparing the simulations in boxes 50 x 50 x 50 s.u. and 
50 X 50 X 100 S.U., one can conclude that interaction with the 
boundary is probably not the main factor in the conversion, as 
it takes significantly different times in the boxes of the same 
width, and similar times in boxes of different widths. 

Critical curvature is a more plausible cause for conversion. 
On qualitative level, larger means smaller initial curvature 
and, by Eq. ([3]), slower evolution. This is indeed what hap- 
pens: the helix with = 100 expands more slowly than 
those with 



50. See Fig. 15 a). Thus it will take longer 
for the filament curvature to reach any particular value, e.g. its 
critical value. 

To quantify this argument, consider an ideal helix whose 
curvature k depends on the radius of its projection p as 



(8) 



p2 + {Lzl2^Y 



where is the height of the helix making one full turn. Fig- 
ure [TSjb) shows how the curvature k of an ideal helix changes 
with p, for the selected values of L^. The important feature 
is that the curvature of a helix increases with p up to a certain 
maximum which depends on . The larger L^, the slower the 
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FIG. 13. (Color online) Same as in Fig.[T2]except with a box of size 
50 X 50 X 100 s.u. 
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FIG. 15. (Color online) Curvature evolution, (a) The mean radius, 
p, of the (x, y) projection of the filament, against the period number, 
Nt. The vertical lines separate the slow region (to the left of the line) 
from the fast region (to the right), (b) Dependence of the curvature, 
on the radius of the projection, p, using Eq. ^ for an ideal helix. 
The horizontal line is the critical curvature, A^crit- 



from an ideal helix. 

We have explicitly calculated the curvature of the filaments 
using numerical differentiation with Tikhonov regularization. 
The details of the procedure can be found in Appendix A. 
The results of this analysis are presented in Fig. [T6ta) and 



Fig.[T7[a), with the local curvature shown against z coordinate 
on the filament, for selected values of the scroll period num- 
bers Nt. We see that the curvature grows very non-uniformly, 
both along the filament and in time. 

To detect the conversion visually, so as to be more precise 
about time and location of the conversion, we have analyzed 
the fine structure of the scroll cores, using a modified version 
of the instant phase singularities: = —1.04, = —0.656. 
The idea is that the point (u^^v^) is within the small loop 
of the phase trajectory shown in Fig. [TJc), and existence of 
such a loop is the signature of slow scrolls. Using this defi- 
nition, slow scrolls are characterized by a double singularity 
line, whereas the fast scrolls have a single singularity line. 



(h) 



See Fig.[T6[b-d) and Fig.[T7[b-d). One can see that the con- 
version is indicated by the departing of the secondary singu- 
larity from the main one following the period in which the 
filament curvature (obtained for the standard, "robust" singu- 
larity = = 0) has exceeded threshold on a substantial 
continuous interval. We see this as a confirmation that cur- 
vature of the filament is the likely cause of the conversion in 
these simulations. 



FIG. 14. (Color online) Same as in Fig.[T2]except with a box of size 
70 X 70 X 50 s.u. 



C. Dynamics of alternative ring vortices 



growth of Hi with p and the smaller the maximum. For com- 
parison, we also show on Fig.[T5jb) the curvature A^crit, which 
causes immediate conversion. This is known from the 2D sim- 
ulation shown in Fig.jSJc) where £ = 0.06 causes immediate 
conversion. Again using the formal equivalence between elec- 
trophoresis and curvature (Sec. |IIBl ), we have A^crit = 0.06. 

These graphs imply that for an ideal vortex, the conver- 
sion due to curvature will happen later for larger L^, which 
is in qualitative agreement with the observations shown in 
Fig.[T5ja). However, there is no quantitative agreement. This 
is likely because the shape of the filaments very soon deviates 



We have also considered the dynamics of alternative scroll 
rings with filament tension of opposite sign at (3 = 0.71. One 
expects that, at least initially, a ring with positive tension will 
drift and contract, while a ring with negative filament tension 
will drift and expand. We are interested in the long-term dy- 
namics of such rings. 

A quarter of a scroll ring with positive filament tension is 
initiated in a cubical box of size 50 x 50 x 50 s.u. Two cases 
have been considered. In one the box has Neumann boundary 
conditions on all sides (not shown) and in the other the box 
has Neumann boundary conditions on four sides and periodic 
boundary conditions on the top and bottom (case shown in 
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FIG. 16. (Color online) Curvature and conversion in the 70 x 70 x 50 s.u. domain corresponding to the simulation shown in Fi^ 
curvature against z coordinate on the filament, for selected values of the scroll period numbers Nt- (b,c,d): the double- singularity 
of the scrolls (b) immediately before, (c) in the very beginning, and (d) after completion of the conversion. 
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FIG. 17. (Color online) Same as in Fig. [16] except for the 50 x 50 x 50 s.u. simulation shown in Fig. [13] 
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FIG. 18. (Color online) Evolution of a fast scroll ring with positive 
tension at /3 = 0.71. One quarter of the ring is simulated in a do- 
main with Neumann boundary conditions on four sides and periodic 
boundary conditions one the top and bottom. Visualizations and time 
series are as in previous figures. 



Fig. 18). The initial radius of the ring is ^ 25 s.u. and meets 



two adjacent sides at right angles as seen in Fig.[T8]^a). 

The scroll drifted upwards, remaining essentially flat. In 
the case with Neumann boundary conditions on all sides, 
the ring approaches the upper boundary where it shrinks and 
eventually disappears in a top corner of the box. With peri- 
odic boundary condition, as seen in Fig. [18] instead of col- 
lapsing, the ring drifts upwards continuously and reaches an 
asymptotic constant radius. In a way, this perpetual move- 
ment is similar to what was observed in 2D simulations 
of electrophoretic drift with periodical boundary conditions, 
Fig.[7];d-f). 



The existence of the stable vortex rings in excitable media 
was first reported by Winfree"^^. The upward component of 
the ring's velocity at the chosen set of model parameters is in 
agreement with the asymptotic theory which gives C3 < for 
(3 = 0.71. More importantly, it seems that the positive fila- 
ment tension, 62 > 0, on its own may be not enough for a ring 
to collapse. It seems that, in the absence of other perturba- 
tions, the positive filament tension shrinks a ring vortex just 
to a stable minimum radius not equal to zero, while the inter- 
nal interaction of parts of the vortex prevents it from complete 
collapse. Note that in simulations of vortex rings with positive 
filament tension, with both Neumann and periodic boundary 
conditions, there no conversion of the fast vortex into its slow 
counterpart has been observed. 

Now consider the evolution of the slow vortex ring with 
negative filament tension shown in Fig. [19] The model pa- 



rameters are the same as in the previous simulation although 
here we show the case with Neumann boundary conditions 
on all sides of the domain. In accordance with the negative 



filament tension, the ring initially expands [Fig. 19 a-b)] and 
becomes non-planar. After five rotations the vortex converts 
into its fast counterpart (with positive tension), as seen in the 



action potential recordings in Fig. [T9];g). From this point the 
vortex propagates to the top boundary where it contracts and 
disappears in the top corner, Fig.[T9l^c-f). 



In case of periodic boundary conditions (not shown) the 
scenario is: expansion due to negative tension followed by 
conversion to the fast vortex with positive filament tension. 
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FIG. 19. (Color online) Evolution of a slow scroll ring with negative 
tension at /3 = 0.71. The simulation is initiate with a quarter of 
the ring a domain with Neumann boundary conditions on all sides. 
Visualizations and time series are as in previous figures. 



followed by contraction to a stable ring with a nonzero ra- 
dius, similar to Fig. The same type of scenario was first 
reported by Sutcliffe and Winfree^^ who explained it by the 
twist of the filament. Here we have demonstrated that the 



switching from expansion to contraction could be due to sim- 
ple switching from a vortex with negative filament tension to 
its counterpart with positive filament tension. 



D. Alternative vortices with negative filament tension 
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FIG. 20. (Color online) Evolution when alternative scrolls both have 
negative filament tension, /3 = 0.73. Simulations start with a fast he- 
lical scroll. The boundary conditions are Neumann on all boundaries, 
(a-f) (Top row) isosurfaces of the i^-field; (second row) filaments 
only, (g) Action potentials with blue crosses marking the times at 
which the snapshots are taken, (h) Action potentials for a simula- 
tion with the same initial and boundary conditions, but with a shock 
stimulus applied at the instant marked with a blue dot). 



At 13 = 0.73, the alternative solutions both have negative 
filament tension, as can be seen Fig. |3] Still, only the slower 
vortex displays DADs in its action potentials. As we know 
from the 2D studies [Fig.|5|, at these model parameters con- 
version is possible both ways, from slow to fast and from fast 
to slow vortex. 

First we consider the evolution of a helical vortex initiated 



from the fast spiral shown in Fig. 20 Simulations are in a box 
of size 50 X 50 X 50 s.u., with Neumann boundary conditions 
on all sides. The filament initially expands due to negative 
tension, but quickly converts to the slow scroll after just two 
rotations. See the DADs appeared in Fig.[20tg). The slower 



vortex continues to expand and then converts back to the fast 
vortex after 14 periods. Since both alternative vortices have 
negative tension, eventually full fibrillation with multiple fil- 



aments develops as is seen Fig. 20 c). There are no further 
spontaneous conversions, and the fast vortex with negative fil- 
ament tension continues to evolve. However, we see a gradual 
reduction in the number of filaments until only a single tiny 
piece persists in the bottom right corner of the box. 

Thus far, we have observed only spontaneous conversion 
between alternative vortices in 3D. It is important for cardio- 
logical applications to verify whether 3D vortices can be con- 
verted by a shock. Having the prolonged period of slow vortex 
evolution in the simulation shown in Fig. [2Q|g), we have re- 
peated the simulation and this time applying a uniform shock, 
Eq. (|7]), to the slow vortex. The action potential time series, 
but not the visualization, are shown in Fig.[20jh). A shock of 
amplitude A = 0.8 is applied just after the third slow vortex 
period. The shock successfully converts the slow vortex into 
its fast counterpart. The rest of the evolution of the fast scroll 



was similar to Fig. 20 a-f): development of the full fibrillation 



followed by gradual reduction of the number of filaments in 
the box to just a single tiny persistent piece in the top right 
comer of the box, but in approximately half of the time taken 
without the shock. 



V. DISCUSSION 

In this study we have investigated dynamics of alternative 
spiral and scroll waves in the FitzHugh-Nagumo model, in 
a parameter region of bistability between distinct alternative 
vortices. Some of the features have been noted already by 
Winfre^who discovered such alternative solutions. Here we 
have extended those studies and related them to the asymp- 
totic theory of spiral and scroll wave dynamics. The most 
important features are: 

1. There is a parameter region in FHN model where 
stable alternative spiral wave solutions, with differ- 
ent rotation frequencies, exist. Faster spiral waves 
have the normal morphology of the "action potential" 
whereas the slower spiral waves show "delayed after- 
depolarization" morphology. 

2. The alternative vortices can convert one into the other 
as a result of a variety of perturbations, such as a homo- 
geneous pulse stimulus or interaction with a boundary. 
There are distinct parametric regions where this transi- 
tion may occur in one direction only, either from slow 
to fast vortex or from fast to slow, as well as the re- 
gions where the transition can occur in both directions. 
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Conversion of a single slow vortex into multiple fast 
vortices is also possible. 

3. The conversion effects of spirals can also be observed 
for scroll waves. In addition, scroll waves can demon- 
strate spontaneous conversion apparently related to the 
filament curvature. 

4. In accordance with the predictions of the asymptotic 
theory, the alternative scroll waves have significantly 
different filament tensions. In particular, the filament 
tension of the faster scroll waves changes sign in the 
parametric interval considered, whereas the slow scroll 
waves have negative filament tension in this interval. 
Therefore, there exists a parameter region where the fast 
vortex with positive filament tension has a slow coun- 
terpart with negative filament tension, and a parameter 
region where both alternative vortices have negative fil- 
ament tension. The scrolls with positive filament ten- 
sion have tendency to contract, or straighten up if the 
filament connects opposite surfaces of the box. The 
vortex rings with positive filament tension might have 
a stable radius not equal to zero, and collapse only due 
to an additional perturbation, e.g. hitting the boundary 
perpendicular to the axis of the ring. 

The scrolls with negative filament tension have the ten- 
dency to lengthen and multiply, which can lead to "scroll wave 
turbulence", phenomenologically similar to certain stages of 
cardiac fibrillation. However, interaction of different scroll 
filaments with each other and/or with the boundary may lead 
to the stabilization of the scrolls notwithstanding the effects of 
the filament tension, as was observed previously (e.g. Ref.[28l 
p. 139) and in simulations presented here. Therefore, in 
three dimensional experiments and simulations the stabilized 
scrolls with negative filament tension may behave identically 
to the scrolls with the positive filament tension, and the only 
way to distinguish between them is to compute/measure their 
response functions. 

The conversion processes summarized in points [2] and |3] 
above are essentially threshold effects and so are not asymp- 
totic in nature. However, there are examples where (non- 
asymptotic) analytical approaches have been successful in 
describing threshold phenomena,^. Following the idea of 
Ref. 51, it may be possible to describe conversion of spiral 
waves in response to external stimuli using center- stable space 
of the unstable spiral wave solutions, which presumably sep- 
arates (in the functional space) the stable alternative spirals 
observed in simulations. 

Assuming that the above features are present in other mod- 
els, including more physiologically realistic ones, one may 
conjecture the following scenarios, which may be relevant to 
cardiac electrophysiology. 

Shock-induced conversion of fibrillation to tachycardia. 
Suppose the cardiac tissue fibrillates due to "scroll turbu- 
lence" mechanism, underlied by slow scrolls. Such slow 
scrolls are likely to convert to fast scrolls either due to cur- 
vature or to interaction with each other or with boundaries. 
However this may take a long time. An electric shock which 



is too weak to instantly defibrillate, may still be enough to 
convert from slow to fast scrolls. If the fast scrolls have posi- 
tive filament tension, this may lead to "delayed defibrillation", 
when the fibrillation stops via collapse of all scrolls, but only 
many cycles after the shock, or to "tachycardia" when the fast 
scrolls would stabilize either by connecting opposite surfaces 
of the tissue (say transmural filaments) or by attaching to lo- 
calized inhomogeneities or anatomical features. 

If the corresponding fast scrolls have negative filament ten- 
sion, they might stabilize due to interaction with each other 
and with the boundary, leading to tachycardia indistinguish- 
able from the one with positive filament tension and without 
pinning to anatomical obstacles. 

Intermittent fibrillation. If both alternative scrolls have neg- 
ative filament tension and can be converted equally one into 
another, then a relatively weak external shock applied to the 
stabilized fast scroll can convert it back into the slow one, 
which will initiate another, may be prolonged fibrillation pe- 
riod before eventual conversion and stabilization back into 
tachycardia. 

FitzHugh-Nagumo model is admittedly very far from real- 
istic cardiac models. However, delayed after-depolarization 
is indeed known in electr ophysiology, and known to have ar- 
rhythmogenic The presence of DADs will add to 

the period of the vortices, and the tendency of vortices with 
longer periods to have negative filament tension is universal^^' 
and not restricted to the FHN model. The coexistence of vor- 
tices with different periods and different filament tensions is 
key to the effects we have described. Hence we believe that 
the basic phenomenology involved in our results is not re- 
stricted to FHN and may be observed in any excitable media 
with DADs, and so is likely to be relevant to cardiac electro- 
physiology. This capacity of DAD is distinct from its well- 
known role as a mechanism of arrhythmia initiation. 
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